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Abstract 

We describe the viscous evolution of a collinear three vortex structure that corresponds initially to an 
inviscid point vortex fixed equilibrium, with the goal of elucidating some of the main transient dynamical 
features of the flow. Using a multi-Gaussian 'core-growth' type of model, we show that the system 
immediately begins to rotate unsteadily, a mechanism we attribute to a 'viscously induced' instability. 
We then examine in detail the qualitative and quantitative evolution of the system as it evolves toward 
the long-time asymptotic Lamb-Oseen state, showing the sequence of topological bifurcations that occur 
both in a fixed reference frame, and in an appropriately chosen rotating reference frame. The evolution 
of passive particles in this viscously evolving flow is shown and interpreted in relation to these evolving 
streamline patterns. 

1 Introduction 

When point vortex equilibria of the 2D Eulcr equations (inviscid) are used as initial conditions for the corre- 
sponding Navicr-Stokes equations (viscous), typically an interesting and complex dynamical process unfolds 
at short and intermediate time scales, which depends crucially on the details of the initial configuration. 
For long enough times, Gallay & Wayne proved recently that the Lamb-Oseen solution is an asymptotically 
stable attracting solution for all (integrable) initial vorticity fields [T] . While very powerful, this asymptotic 
result does not elucidate the intermediate dynamics that take place in finite time and allow a given initial 
vorticity field to reach the single-peaked Gaussian distribution of the Lamb-Oseen solution. Given the rather 
large (and growing) literature on point vortex equilibria of the Euler equations, (see for example, [3[3l|4]), 
we thought an analysis of how these equilibria evolve under the evolution of the full Navier-Stokes system 
would merit a systematic treatment. Hence, in this paper, we begin an investigation of the viscous evolu- 
tion of a class of initial vorticity fields consisting of the superposition of N Dirac-delta functions or point 
vortices [2J. Our initial configuration, shown in Figure [Tj is a collinear configuration of three point vortices, 
evenly spaced along a line (the x axis), with strengths 2r, —T, and 2r, respectively. Such a configuration, 
for the Euler equations, is known to be an unstable fixed equilibrium, as fleshed out most recently and com- 
prehensively in [5], but earlier in [BJ. We point out that this configuration, because of the strengths chosen 
for each of the point vortices, is not what is commonly referred to as the 'tripole' state [71 [HI ESI EH] m which 
the vortex strengths sum to zero. Our focus in this paper will be the dynamics that takes place at the short 
and intermediate time scales, using this initial state in the Navier-Stokes equations, before the long-time 
asymptotic Lamb-Oseen solution dominates. This includes the dynamics of the surrounding passive field 
and the corresponding background time-dependent streamline pattern in an appropriately chosen reference 
frame which we argue is very helpful as a diagnostic tool to interpret the resulting flowfield. 

2 Problem Setting 

Consider an incompressible fluid in an unbounded two-dimensional (2D) domain R 2 . The fluid motion is 
governed by Navier-Stokes equations, written in terms of the vorticity field w(x, t), a scalar- valued function 
of position x and time t, as follows 

— = -u • Vw + isAu . (1) 
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Figure 1: Fixed point- vortex equilibrium: three collinear and equally-spaced point vortices with the outer vortices of strength 
2r and the middle vortex of strength — T. 



The kinematic viscosity v is assumed to be constant. The fluid velocity u(x, t) is a vector-valued function of 
x and t. Both x and u are expressed in an inertial frame {ei}i=i,2,3 where (ei, e2) span the plane of motion, 
that is to say, one has x = x ei + y e 2 and u = u x ei + u y e 2 or, equivalently, x = (x, y) and u = (u x , u y ). 
By definition, the vorticity vector u> — V x u is always perpendicular to the plane of motion and can thus 
be expressed as iv = lu e,3 . The velocity u and vorticity uj are related via the 2D Biot-Savart law 

u(x, t) = ~ / [ X - X ^ 9 w(x, t) dx , (2) 
2ti" Jg? ||x - x||- 

where x is an integration variable and x = (— y, x). Note that for 2D flows, the stretching term u> ■ Vu 
is identically zero thus does not appear in ([I]) while the continuity equation div(u) = is trivially satisfied 
when expressed in terms of vorticity. 

The solution of the system of equations ([I]) and ^ depends, of course, on the choice of initial 
conditions cj(x, 0). One solution of particular interest in this work is the well-known Lamb-Oseen solution 
corresponding to a Dirac-delta initial condition w(x, 0) = r<5(x), i.e., a point vortex placed at the origin with 
circulation or strength F (more generally, the circulation F around any closed curve C in the fluid domain is 
defined as Y = § c u • ds = J A ujda and can be thought of as the flux of vorticity through the area A enclosed 
by the curve C). Traditionally, the problem can be expressed compactly in complex notation with position 
variable z, z = x + iy and i = \J — 1. The Lamb-Oseen solution is given by (see, for example, |lll, [2"] for more 
details) 

w(M) = exp ( -^r-} • (3) 



(4) 
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The corresponding velocity field u is given by, 
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m v = 7^—-- 
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1 — exp 



where the dot notation () = d()/dt refers to the time derivative and the notation z* = x — iy refers to the 
complex conjugate. According to Q, the evolution of a vorticity field that is initially concentrated at the 
origin is such that the vorticity diffuses axisymmetrically as a Gaussian distribution. The spreading of the 
vorticity concentration can be quantified by the vortex core or support defined as p — v 1 4vt — \/4t, where 
r = vt. 

Despite the explicit, exact, and simple nature of the solution (13]), for more complicated initial data, 
explicit solutions of and ^ are not analytically available for a general initial vorticity field. We are 
particularly interested in the viscous evolution of a class of initial vorticity fields w(z, 0) = -, r Q <5(z — z Q ) 
consisting of the superposition of N Dirac-delta functions or point vortices. Note that Gallay & Wayne proved 
recently that the Lamb-Oseen solution is an asymptotically stable attracting solution for all (integrable) 
initial vorticity fields, see [T], of which w(z, 0) = 53 Q =i -T Q (5(z — z a ) is a special case. It is the dynamics that 
unfolds as the system evolves towards this final state that we are interested in. 

The dynamics of N point vortices in the plane is extensively analyzed in the context of the inviscid 
fluid model (y — 0), see for example j2] and references therein. The vorticity field remains then concentrated 
for all times at N points whose position z a (t) (a = 1, . . . , N) is dictated by the local fluid velocity induced 
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by the presence of the other vortices. The fluid velocity at an arbitrary point z in the plane that does not 
coincide with a point vortex is obtained from the Biot-Savart law ([2| which takes the form 



A i r 



^— ' 2-7TZ Z — Z c 



(5) 



whereas the velocity at a point vortex zp is given by subtracting the effect of that point vortex from d5]), 
and replacing z with zp, namely, 

N 

**B = Y7T — • (6) 

£-? 2m zp-z a K ' 

The 2N first-order ordinary differential equations ^ dictating the inviscid evolution of N point vortices are 
known to exhibit regular, including fixed and moving equilibria, as well as chaotic dynamics depending on 
the number of vortices, their strengths and initial positions. The literature on this general topic is large, and 
we refer simply to the influential 1983 review article of Aref [T^], along with the monographs of Saffman [TT] 
(especially Chapter 7) and Newton [5] for an immediate entry into the literature. We also mention the 2008 
IUTAM Symposium '150 Years of Vortex Dynamics' held at the Technical University of Denmark in which 
the lively state-of-the-art developments were reported [15] . 

In order to highlight the way in which the presence of viscosity affects the inviscid point vortex 
dynamics, we focus on studying the viscous evolution of N point vortices whose location at time t = 
correspond to a fixed equilibrium of the inviscid point vortex model ([6]). For concreteness, we consider the 
case of N = 3 collinear and equally-spaced point vortices as shown in Figure [l] Let ro denote the distance 
between two adjacent vortices and let 2r be the circulation of the left and right vortices and while — T be 
that of the center vortex. Also, let zl, zq and zr denote the positions of centers of the left, center and right 
vortices, respectively. One can readily verify that this configuration constitutes a fixed equilibrium of the 
inviscid point vortex model Q. 

It is convenient to non-dimensionalize the problem using the length scale L = 2r§ and the time scale 
T dictated by T, namely, T — L 2 /V . To this end, the non-dimensional parameters are given by 

f o=2> f = l, v=f, ( 7 ) 
and the non-dimensional variables are of the form 

Z _ U bJ 

Z= L' U= L]f' ^Ijf- (8) 

For simplicity, we drop the ~ notation with the understanding that all parameters and variables are non- 
dimensional hereafter. The Navier-Stokes equation ([7J) can be rewritten in dimensionless form 

dbJ „ 1 . 

— = -u-Vw + — - Aw, 9 
ot Re 

where Re is the dimensionless Reynolds number, here defined as Re = T/u. Some words of caution are in 
order here, as we will be comparing numerical simulations of the Navier-Stokes equations with our model, 
and to do so requires that one is able to compare the direct numerical simulation (DNS) Reynolds number 
with the 'model' Reynolds number. For this, it is better to think of the Reynolds number as the ratio of 
inertia! effects — u • Vw over diffusive effects Aw. In some sense, one can think of the term — u • Vw as being 
primarily responsible for the rotation we will discuss, while the term Aw not only triggers the rotation, but 
diffuses the cores of the vortices. For any DNS, this creates an 'effective' Reynolds number which depends 
not only on the choice T/u, but also on numerical discretization since it affects the 'rotation' and 'diffusion'. 
The 'model' Reynolds number will be discussed in more detail in the upcoming sections. 

By way of motivation, we first present a numerical solution of the system of equations (|9| and ^ 
subject to the initial condition w(z, 0) = 2TS(z — ro) — T(5(z) + 2r<5(z + ro). We use the numerical algorithm 
devised in 14J that utilizes a second-order finite difference method with a multi-domain non-reflecting bound- 
ary condition to emulate the infinite fluid domain. This is a mesh-based method which poses a problem in 
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(a) t = (b) t = 2.8 (c) t = 43.7 (d) t = 47 .4 




(c) t = (f) t = 2.8 (g) t = 43.7 (h) t = 47.4 

Figure 2: Vorticity contours (top row) and streamlines (bottom row) of Navier-Stokes simulation for Re = 1000 at t = 
0,2.8,43.7 and 47.4. The vortex configuration rotates unsteadily for t > 0. The center vortex stretches and diffuses out first, 
then the outer two vortices begin to merge. Eventually the vortex configuration approach a single Gaussian vortex. 



handling the Dirac-delta initial conditions because they are not well-posed for discretization on a standard 
Euclidean mesh. To overcome this problem, we consider the initial conditions of the vorticity field as a super- 
position of three slightly diffused Gaussian peaks. In all the simulations presented here, we diffuse the initial 
Dirac-delta vorticity field by e such that e/Re — 2 x 10~ 4 , e.g., e = 0.2 when Re = 1000. This, of course, 
introduces a slight mismatch in the initial conditions used for the numerical simulation with those used in the 
model, an error we cannot completely eliminate, but should be kept in mind when comparing the simulation 
with the model. We compute the time evolution of the vorticity field in the window [—1.5, 1.5] x [—1.5, 1.5] 
while the non-reflecting boundary conditions are imposed using the multi-domain technique with 10 nested 
domains, the largest of which is 2 10 times the size of our result window. The spatial and time steps are set 
to Ax = Ay = 0.01, At = 0.02. 

Figure [2] depicts the time evolution of the vorticity contours (top row) and streamlines (bottom row) 
of Navier-Stokes solution for Re — 1000 at four time instants: t = 0, 2.8, 43.7 and 47.4. One notices that 
the vortex configuration begins to rotate unsteadily for t > 0; we refer to this motion as viscosity-induced 
rotation. One also notices that the center vortex stretches and diffuses out first, then the outer two vortices 
begin to merge. Eventually the vortex configuration approaches a single Gaussian vortex. 

3 The multi-Gaussian model 

In this section, we use a simple, analytically-tractable model to describe the dynamic evolution of N point 
vortices for non-zero (but small) viscosity (y ^ 0). The model assumes that the vorticity of each initial 
point vortex spreads axisymmetrically as an isolated Lamb-Oseen vortex, thus modeling the diffusion term 
vAlj in ([lj, while its center moves according to the local velocity induced by the presence of the other 
(diffusing) vortices, thus accounting for the convection term u • Vw in |l]). It is worth noting here the recent 
work of Gallay [H] who analyses the inviscid limit v — > of the 2D Navier-Stokes evolution of Dirac-delta 
initial conditions and proves, under certain assumptions, that the solution of the Navier-Stokes equation 
converges, as v — > 0, to a superposition of Lamb-Oseen vortices. In this work, we show that, for small yet 
finite is, the multi-Gaussian model is able to capture qualitatively, though not quantitatively, some of the 
main features of the Navier-Stokes solution. Generally speaking, this class of models has been most highly 
developed in the numerical literature (see for example [THIE] and subsequent analysis in [TH]) and is referred 
to as a 'core-growth' class of models. One can trace the 'splitting' idea of the advection and the diffusion 
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(a) t = (b) t = 25 (c) t = 60 (d) t = 300 




(c) t = (f ) t = 25 (g) t = 60 (h) t = 300 



Figure 3: Vorticity contours (top row) and streamlines (bottom row) of multi-Gaussian model for v = 1/1000 and T = 1 at 
four instants t = 0, 25, 60 and 300. Similar to Navier-Stokes simulation, the vortex configuration rotates unsteadily for t > 0, 
the center vortex stretches and diffuses out first, then the outer vortices merge, eventually the vortex configuration approaches 
a single Gaussian vortex. 



terms of the 2D Navier-Stokes equations, on which the core growth model is based, at least back to Chorin's 
influential paper [19], also used by Milinazzo and Saffman [20]. In these papers, the diffusion was handled by 
a random walk approach. Core growth models based on time-dependent solutions of the heat-equation were 
developed and used mostly by the numerical/computational vortex dynamics community, and are discussed 
and developed explicitly in [2U 123 • In the context of numerical simulations, focused studies can be 
found in the works of Barba and Leonard [25] , and used in specific models in [57J [21] • We mention, 
of course, also the works [30l EU Q] and the 2009 Ph.D. thesis of Uminsky [32] and follow-up work [33] 
which develops an eigenfunction expansion method based on the form of the heat-kernel. Additionally, we 
mention the body of work generated by Dritschel and co-workers, of which [34, 35 would be two relevant 
examples, whose aim is to elucidate via Lagrangian type numerical simulations, the host of complex processes 
associated with mixing and dynamics in viscously evolving two-dimensional flows. 

The model assumes that the vorticity field at all times is a superposition of multiple Lamb-Oseen 
vortices 

= f: -fe exP (±=#^y do) 



The associated velocity field is computed by substituting ( 10 1 into |2]) 



\ -> 1 a 

2m(z - z a ) 



1 — exp 



2 \ 1 



(11) 



The velocity at the center of the [3 th vortex is given by subtracting the effect of that vortex and replacing 
z by zp in (fTTJ) 
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E 



27ri(z /3 - Z a ) 



1 — exp 



(12) 



The system of equations in ( 10 1, ( 11 1, and ( 12 1 is referred to as the multi- Gaussian model. We emphasize 
that we include in this model the equation (11) for the evolution of passive tracers in the 'field' which is 
transported under the dynamics generated by (10) and (12). This will be discussed more thoroughly in 
Section [5] and is relevant for comparisons of panels (a)-(d) of Figure [2] with (a)-(d) of Figure [3] 
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Figure 4: Rotation rate 8 and rotation angle 9 as functions of time r = vt of multi-Gaussian model for T = l,ro = 0.5, v 
10" 3 . 



According to (10), the vorticity field associated with the initial three- vortex configuration shown in 
Figure [T] is given by 



u(z,t) = 



2Y exp 



Avt 



Texp 
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(13) 



The location of the centers of the vortices z^, z^ and z^ is obtained by solving the set of six first-order, 



ordinary differential equations in (12). From symmetry, one can readily verify that i* c = and that the 



centers of the vortices remain collinear and equally-spaced with constant distances for all time. The vorticity 
contours of (10) are then plotted in Figure [3] (top row). The streamlines associated with the velocity field 
in (12) is shown in Figure [3] (bottom row). Similarly to the Navier-Stokes solution depicted in Figure[2j the 
dynamic evolution of the multi-Gaussian model is characterized by: (i) an unsteady rotation of the whole 
vortex configuration for t > 0, (ii) a stretching of the middle vortex, and (iii) eventual merging of the outer 
two vortices to form one single-peaked Gaussian of strength 3r as shown in Figure [3] However, here some 
care is in order, as clearly Figures [2](b)-(d) (DNS) and Figures §»-(d) show some important differences. 
Not only are the timescales different, but Figure [2jb) shows a convective 'wrappping' and 'stretching' of 
the middle vortex around the outer two before the diffusive effects kick in, whereas Figure [3^b) shows the 
stretching, but not the wrapping. Here it is important to remember that the passively advected field, as 
shown in Figure [l3j is an important part of the model, and this field does show some of the same nonlinear 
wrapping features that appear in the DNS Figure [2jb). One could say, in some respects, that the outer two 
vortices, being twice the strength of the inner one, are the primary 'drivers' of the flowfield, which is perhaps 
why Figure [2je)-(h) match relatively well with Figure[3](e)-(h). The 'passively advected' inner vortex shown 
in Figure [2jb) is better reflected in the passive particle field shown in Figure 13 and discussed at length in 
Section [5] In turn, because the passively advected field in our model is not affecting the vorticity evolution, 
whereas in the DNS it is, this helps explain why the timescales associated with the two are different. The 
model is not an exact solution of the Navier-Stokes equations, and this appears to be the main physical 
manifestation of this fact. 



The unsteady rotation rate of the vortex structure is obtained analytically as follows. From (12), the 
velocity of one of the outer vortices, say the right vortex, takes the form 
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(14) 



Now, by symmetry one has z^ = —roe , zp = and zr = 
traced by the vortex centers and the x axis, and z# = ir§Be l 
simplifying, that the rotation rate 9 is given by 



roe , where 9 is the angle between the line 



One gets, upon substituting into (14) and 
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(15) 



In Figure [4] is a depiction of 9 versus t — vt which shows that the rotation rate starts from zero, reaches a 



maximum value 9 rnax at an intermediate time r„ 



vt„ 



and eventually decays to zero as vt —> 00. The 
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(a) t = 



(b) r = tl fa 0.016 



(c) t = 0.05 



(d) t = t 2 as 0.086 



Figure 5: Evolution of the streamlines of the multi-Gaussian model. The separatrices are depicted in thick lines with arrows 
showing the direction of the flow. Instantaneous hyperbolic points are at intersections of separatrices while elliptic points are 
represented by circles. 
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Fi glire 6; Instantaneous stagnation points: (a) a pair of hyperbolic stagnation point located at (0, i??/) for r = 1,7*0 — 0.5 
and v = 10" 3 . This pair collides at (0, 0) at bifurcation time t\ ~ 0.016. (b) a pair of elliptic stagnation points located at 
(±£y,0). This pair collides with the now hyperbolic origin at bifurcation time T2 ~ 0.086. 
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are given by 
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The orientation angle 9 can be readily obtained by integrating ( 15 ) in time 
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(16) 



(17) 



where the exponential integral is defined as Ei(x) = — Jl^ exp(t)/t dt in the sense of principle value, which 
can be evaluated numerically to machine accuracy. 

It is convenient for analyzing the flow to express the fluid velocity field z in a frame co-rotating with 
the vortex configuration at the time-dependent rotation 9. Let £ = (C, J?) denote position of a point in the 
plane expressed in the rotating frame. The point transformation from the rotating to the inertial frame is 
given by 

r cos 9 — sin f 
sin 9 cos ( 



z = R£, R = 



(18) 



The fluid velocity transforms as z = Rv where v = (vq,v v ) is the velocity field expressed in the rotating 
frame. In component form, equations (11) transform as 
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(19) 
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(20) 



where we used r = vt. 

The instantaneous stagnation points of the velocity field (obtained by setting the right-hand side of 



(191, (20) to zero) reveal important information about the instantaneous streamlines of the fluid velocity 



field. From symmetry of the velocity field, the instantaneous stagnation points must lie on the £ and r\ axes. 
One finds a total of five fixed points: one initially elliptic point at the origin, a pair of initially hyperbolic 
points at (0, ±77/) and a pair of initially elliptic points at (±£/,0). The hyperbolic and elliptic character 
of these stagnation points is obtained by linearizing (l9| ,(20) about the instantaneous stagnation points 
and computing the eigenvalues of the linearized system. A pair of eigenvalues ±A is associated with each 
stagnation point. One has a hyperbolic point if A is real and an elliptic point if A is pure imaginary. 

At time t = 0, the streamlines are those of the inviscid equilibrium, with a separatrix linking the 
two hyperbolic stagnations points on the 77 axis, as shown in Figure [5j Initially, the separatrix divides the 
fluid domain into four regions: three regions, one around each point vortex or elliptic point, and a fourth 
region bounded by the separatrix and the bound at infinity and void of point vortices. As time evolves, 
the location of the instantaneous stagnation points change as shown in Figure [5] and the separatrix evolves 
accordingly. Note that the time-dependent separatrix does not constitute barriers to fluid motion and fluid 
particles typically move across this separatrix as time evolves as discussed in more details in Section [5] 
Figure [6] shows the coordinates of the stagnation points ±77/ and ±C/ as functions of time. The pair of 
initially hyperbolic points (0, ±77/) start from (0, ±ro/y/S) then collide together with the elliptic point at the 
origin in finite time t\ ~ 0.016 to transform the origin into a hyperbolic point. This collision of instantaneous 
stagnation points is accompanied by a change in the streamline topology where the region around the center 
vortex disappears, see Figure [5j Time t\ is referred to as the first bifurcation time. (Note that the first 
bifurcation time n does not correspond to when the cores of the three Gaussian vortices touch for the first 
time which takes place at r = ro/16 = 0.015625, nor does it correspond to when the cores of the two outer 
Gaussian vortices touch r = Tq/4 = 0.0625. Indeed, the definition of core size of a Gaussian function is 
rather ad hoc and bears little relevance to the dynamics of the multi-Gaussian model.) Similarly, (±£/,0) 
starts from (±r ,0) and collides at the now hyperbolic point at the origin at time t 2 rs 0.086. Time t 2 is 
referred to as the second bifurcation time. For t > r 2 , one has one single elliptic point at the origin as 
expected from the asymptotic Lamb-Oseen solution. 



4 Comparison to Navier- Stokes 



The residual a of the model is computed by substituting the solution of ( 10 ) and ( 11 ) into the Navier-Stokcs 
equation If the solution of the model is also an exact solution of the Navier-Stokes equation for a 
given set of initial conditions, the residual a is identically 0. In general, a is not zero (see discussions of 
this in [21]) and it can be viewed as an indication of the inaccuracy of the multi-Gaussian model. The L 2 
norm of residual is plotted as a function of time r = vt in Figure [7]j a) for the collinear vortex configuration 
considered here. Figure [7ja) shows that as r increases, the L 2 norm of a tends to zero, indicating that 
the multi-Gaussian model agrees with the Navier-Stokes solution for r large. From the result of Gallay & 
Wayne pQ and since the total circulation of the initial vorticity field is 3r 7^ 0, we know as t — > 00 the 
Navier-Stokes solution approaches a single Gaussian vorticity distribution w^, = {2>T/Airvt) exp (— ||z|| 2 /4i/i) 
centered at the origin with circulation 3r. We compute the difference between the multi-Gaussian model 
and this asymptotic solution uj^ . Figure [7](b) and (c) show the L 2 norm of the difference in both vorticity 
and velocity, respectively. These plots confirm that the multi-Gaussian model approaches the asymptotic 
Lamb-Oseen solution for large time but at the intermediate times, the multi-Gaussian model exhibits richer 
dynamics than the asymptotic Lamb-Oseen solution. While the dynamics of the multi-Gaussian model at 
these intermediate time scales does not faithfully track the Navier-Stokes solution (as seen from Figure (7Ta)), 
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Figure 7: Multi-Guassian model: (a) L2 norm of residual a versus r = ut for T = 1 and v = 10 -3 ; (b) L 2 norm of the 
difference in the vorticity field of the multi-Gaussian model and the single peaked Lamb-Oseen vortex with circulation 3F; and 
(c) shows the difference in velocity field. Clearly, for long time, the model approaches the single peaked Gaussian but in short 
time, the multi-Gaussian, while not numerically accurate in comparison to the Navier-Stokes model as indicated in (a), its 
dynamics is richer than the single Gaussian as indicated in (b) and (c). 
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(a) Navier-Stokes simulation (b) Multi-Gaussian model 



Figure 8: Comparison of rotation angle 9 between (a) the Navier-Stokes simulations and (b) the multi-Gaussian model. 
Navier-Stokes simulations are conducted with the same initial vorticity field for Reynolds numbers Re = 1 000, 2000, 3000 and 
4000. The results of the multi-Gaussian model are obtained for T = 1 and v = 1/100, 1/200, 1/300, and 1/400. The trend of 
both models is qualitatively similar. 



it does capture more details than the asymptotic Lamb-Oseen vortex and its evolution seems to exhibit the 
main qualitative features of the Navier-Stokes model as argued next. 

It is evident from Figures [2] and [3] that both the Navier-Stokes equations and the multi-Gaussian model 
exhibit a viscosity-induced rotation as t > 0. In Figure [8] we compare the qualitative trends of rotation angle 
9 obtained from the numerical solution to the Navier-Stokes equation and the analytical solution of multi- 
Gaussian model. In the Navier-Stokes solution, the rotation angle 6 is obtained by computing the angle 
between the line traced by the vorticity peaks (see Figure [2]) and the x axis while it is given by (17 1 in 
the multi-Gaussian model. Clearly, both the Navier-Stokes solution and the model, although quantitatively 
distinct, exhibit similar qualitative trends in that the rotation angle 9 is smaller when Re increases (in 
Navier-Stokes) or equivalently when Y/v increases (in the model). As cautioned earlier about comparing 
DNS Reynolds numbers with the 'model' Reynolds number, if both are thought of strictly as T/u, we can 
only claim qualitative overlap with the model and DNS. To obtain more quantitative overlap would require 
more effort on our part to obtain an accurate Lagrangian based DNS to get a more detailed handle on 
the effective 'numerical' Reynolds number, along with a modified model system that does more to couple 
rotational effects with diffusive effects, neither of which are the immediate goals of the current work. 

To quantify the difference between the Navier-Stokes solution and the multi-Gaussian model, we focus 
on comparing the first bifurcation time T\ in the Navier-Stokes simulation for different Re to the first bifur- 
cation time in the model. The result is plotted in Log-Log scale in Figure [9] for Re = 1000, 2000, 3000, 4000 
and 5000 (plotted in squares) . The dashed line is the best fitted straight line using the least square distance 
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Figure 9: Comparison of the times of the first bifurcation given by the Navier-Stokes simulations and the multi-Gaussian 
model in Log-Log plot. Simulation results are plotted as squares for Re = 100, 500, 1000, 2000, 3000, 4000, 5000 and the dashed 
line is a fitted linear line obtained by least square distance rule. The solid line is the result from Multi-Gaussian model. 



rule. The fitted line can be expressed as ln(ri) = ln(Re) — 5.356, which means in linear scale, the fitting 
is T\ — 0.00472i?e. The simulation results are compared to the first bifurcation time t\ = 0.0l6T/i/ as 
predicted by the multi-Gaussian model. While the first bifurcation in the simulations happens earlier than 
in the model, both the Navier-Stokes simulations and model indicate that T\ is linearly dependent on Re (or 



5 Evolution of vorticity in the multi-Gaussian model 

We now use the multi-Gaussian model to analyze the fluid velocity and vorticity fields at intermediate time 
scales before the asymptotic state of a single Lamb-Oseen vortex dominates. The goal of this analysis is 
to understand the intermediate mechanisms that lead the initial point vortex configuration to reach the 
asymptotic Lamb-Oseen vortex. 

As time evolves, the vorticity field, initially concentrated at = and zl,r — Tl/2, begins to spread 
spatially inducing a velocity field similar to that of a Rankine vortex with time-dependent core. By way of 
background, the reader is reminded that the fluid velocity at a point (£, 77) associated with a Rankine vortex 
at the origin with vorticity 3r (corresponding to the total circulation of the collinear vortex structure) is 
perpendicular to the distance r from the origin and its value is given by 




for r < R r 



for r > R c 



(here r 2 = ( 2 + r? 2 ) . (21) 



The value R cr is referred to as the core of the Rankine vortex. For r < R cr , the fluid velocity corresponds 
to a rigid rotation while for r > R cr , the velocity field decays proportionally to the inverse of the distance 
r. As time evolves, the velocity field induced by the viscously evolving collinear vortex structure becomes 
analogous to that of a Rankine vortex with vorticity 3r and time-dependent core, as seen in Figures [TO} This 
analogy is especially evident in Figures |10( b) where we superimpose the velocity field of the Rankine vortex 
on that induced by the viscously evolving collinear vortex structure at three different instances. Close to the 
origin, the velocity field of the collinear vortex structure looks like a rigid rotation, and the rotation rate is 



given by 6 in (15). Since the rotation rate 9 is unsteady, the core size R cr of the Rankine vo rtex , obtained 
by equating 3T/2irR 2 r = 6, is time dependent and it increases with time t as shown in Figures 10 'a). As the 
distance from the origin increases, the velocity field of the collinear vortex structure decays analogously to 
the inverse decay with vorticity 3F. 

Motivated by this analogy with the Rankine vortex, we examine the time evolution of the relative 
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Fi gure 10; The velocity field induced by the collinear vortex structure becomes analogous to that of a Rankine vortex. 
In particular, the component of velocity u„ along the £-axis is depicted. In (a), we show the velocity profiles for vt = 
0.05, 0.1, 0.15, 0.2, 0.3 and 0.5. The maximum velocity decreases as vt increases. When vt is small, e.g. vt = 0.05, is negative 
close to the origin. This is because the vorticity is still relatively concentrated at the vortex centers. In (b), we superimpose on 
the plots of Vn versus ( (solid lines) the velocity of a Rankine vortex (dashed lines) with vorticity 3r and time-dependent core. 
Clearly, the velocity field is similar to that induced by a rigid rotation close to the origin and it is similar to an inverse decay 
at large distance from the origin. 



velocity field 



i = v - k 1 - 



(22) 



obtained by subtracting a rigid body rotation from the fluid velocity field v expressed in the rotating frame 
(written in ( 19 1, ( 20 1 in component form) . Similarly to the analysis in Section[3j we identify the instantaneous 
stagnation points associated with the relative velocity field ( 22 ) . Immediately as t increases from zero (t > 0) , 



in addition to the elliptic stagnation points located at the origin and (±£/,0) and the hyperbolic points at 
(0, ±?7j), one gets two new pairs of stagnation points appearing from infinity: one elliptic pair located at 
(0, ±17/2) and one hyperbolic pair located at (±£/ 2 , 0). Figure 11 a) shows the values of ±77/2 as functions of 
time. Clearly, ±77/2 start from ±00 and eventually converge to final values ?7/2|t-s.oo = ?'o \J 1 1 /3 ~ 1.9rrj. The 
pair (0, ±77^2) remains elliptic for all time. Figure |TT|(b) shows that the pair (±£y 2 ,0) also starts from ±00 
and reaches ±rrj at 73 ps 0.1904 then the origin at Tg ps 0.2045. Meanwhile, the topology of the streamlines 



of the relative velocity field changes as a result of five distinct bifurcations as depicted in Figure 12 and 
explained next, thus the notation t*,...,t£. 

The first bifurcation in the streamline topology is due to the same mechanism explained in Section [3] 
and takes place at the same time r* = Ti, see Figure [l2^c). The second bifurcation does not coincide in time 
with the second bifurcation identified in Section |3j that is, rj* 7^ T2. It is associated with a change in the 
streamline topology caused by a collapse of the separatrices associated with the hyperbolic pair (±C/2 ; 0) onto 
the separatrices of the now hyperbolic point at the origin, see Figure[l2](e). The third bifurcation occurs at 
when the hyperbolic points at (±C/2, 0) collide with the elliptic points at (±Cf = ± r o, 0), respectively, causing 
them to change to hyperbolic points, see Figure 12 g). After the third bifurcation, one still has two pairs 
(±Cf2j0) and (±£y,0) of stagnation points on the £-axis but with exchanged hyperbolic/elliptic characters. 
The fourth bifurcation takes place at t| due yet to another collapse of the separatrices of the hyperbolic 
point at the origin with the separatrices at the now hyperbolic points at (±ro,0), see Figure 121). The 
fifth bifurcation takes place at Tg ps 0.2045 when the now elliptic pair (±(/2, 0) collides with the hyperbolic 
origin causing it to turn into an elliptic point, see Figure 12 k). This bifurcation sequence turns out to be 



crucial in dictating the time evolution of the vorticity field which we visualize using colored passive tracers 
as commonly done in experimental and computational fluid mechanics (see, for example, [9]). 

We seed the flow at time t = with passive tracers of four different colors as shown in Figure [l3|^ a) to 
distinguish the initial four fluid regions identified in Section [3j namely, the three regions around the vortices 
bounded by the separatrix (seeded with red, blue and green particles, respectively) and the fourth region 
(seeded with yellow particles) bounded by the separatrix and the bound at infinity. We let the passive 
tracers be advected by the fluid velocity field given in (11). Snapshots of the passive tracers at six distinct 
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Fi gure Hi Relative velocity field: two new pairs of stagnation points appear from infinity as time t > 0. (a) 77-component of 
the pair of stagnation points (0,±r;^2)- This pair (0, ±^2) eventually converge to (0, ztro-^/ll/3) as r — > 00, respectively, (b) 
^-component of the pair of stagnation points This pair (±(^2) reach (±ro,0) at bifurcation time 73 ~ 0.1904, and 

collapse at (0, 0) at bifurcation time Tg fti 0.2045. Parameters are r = 1, rg = 0.5 and Re = 1000. 



instants in time are depicted in Figure 13 As time evolves, the location of the stagnation points and the 
associated separatrices change. Due to incompressibility, the particles initially in the region around the 
middle vortex (blue color) 'leak' along the unstable branch of separatrices associated with the instantaneous 
hyperbolic points (0, ±77/2) • At r*, Figure |T3^b) shows that all the particles are squeezed out of the middle 
region. Meanwhile as time progresses, the fluid particles in yellow begin to form lobes that stretch at a 
finite distance away from the initial location of the vortices, see Figure [T3^c). Qualitatively, the passive 
tracers in Figure |13[ c) indicate a vorticity field similar to that obtained from the Navier-Stokes simulation 
in Figure (2jjd) (modulo the rigid rotation of the whole structure) . The formation of these lobes cannot be 
explained based on the analysis of the streamline patterns in Section [3] Indeed, the formation of these lobes 
is initiated when the yellow passive tracers encounter the separatrices associated with the hyperbolic points 
of the relative velocity field ( |22[ ) (±£^2, 0) that appear from infinity and move towards the origin along the 
C-axis (see Figure [lT|b)). The lobes then stretch and rotate around the elliptic points (0, ±77/2) that appear 
from infinity and converge to a finite distance away from the origin (see Figure |TT]( a)). Eventually, the 
passive particles initially placed in the regions around the point vortices, whose detailed evolution is also 
dictated by the sequence of bifurfactions described in Figure [l2j join the large lobes as well and begin to 
stretch and rotate at a finite distance away from the initial vortex configuration, see Figure 13'd)-(f). After 
the last bifurcation in Figure 12 'k), all the passive particles continue to rotate as shown in Figure [l3|[f). We 
emphasize that this interesting dynamics of the passive particles, which in turn indicates the evolution of 
the vorticity field, cannot be explained based solely on the analysis of the streamlines of the fluid velocity 
field of Section [3] In addition, because of the detailed and delicate nature of the full series of topological 
bifurcations that occur, to capture all but the first of these in a DNS would require considerable further 
effort and is beyond the scope of the current manuscript. 



6 Conclusions 

The redistribution (inviscid) and diffusion (viscous) of delta-function initial distributions of vorticity, al- 
though configuration independent for sufficiently long timescales, is highly dependent on the initial positions 
and strengths of the point vortices on short and intermediate timescales. These are typically the timescales in 
which much of the important mixing, transport, and redistribution of vorticity is achieved in many settings. 
Greengard's 1985 paper notwithstanding 21: pointing out that the types of models based on advection and 
core diffusion are not exact solutions of the Navier-Stokes equations, these ideas are exceptionally useful in 
getting a handle on some of the important dynamical mechanisms that occur during the evolution towards 
the ultimate Lamb-Oseen state. In fact, one contribution of the current manuscript is to further quantify 
and understand the limitations of 'core-growth' type models as diagnostic tools for understanding more and 
more complex flows and to point out some of the delicate issues in comparing a DNS with these models. 
Not surprisingly, core-growth type models are also useful as starting points for more sophisticated numerical 
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(i) r = t* w 0.2005 (j) r = 0.202 (k) r = r * « 0.2045 (1) r = 0.23 



Fi gure 12: Evolution of the separatrices of the relative velocity field v — Instantaneous hyperbolic points are at 

intersections of separatrices and elliptic points are represented by circles. The outside separatrices and elliptic points in (b) and 
(c) are plotted out of scale. 



methods which systematically exploit some of the main features [16, 24 (Also, see Blobflow, an open source 
vortex method package developed by Rossi, available at http://www.math.udel.edu/^rossi/BlobFlow as of 
October 2010). 

We summarize here with three main points associated with the viscous evolution of the three-vortex 
collinear state whose initial configuration corresponds to an unstable inviscid fixed equilibrium: 

(i) The presence of viscosity immediately 'triggers' the underlying instability of the equilibrium, caus- 
ing the vortices to rotate unsteadily; 

(ii) In a fixed frame of reference, as the system evolves towards the ultimate Lamb-Oseen solution, 
the streamline patterns associated with the velocity field undergo a clear sequence of topological bifurcations 
which we depict in Figure [14] We show the 'homotopic equivalence' of each of the distinct patterns in 
the panels - the time and quantitative values of the pattern are not depicted, just the sequence of distinct 
patterns that appear during the time sequence; 

(iii) More interestingly, since the velocity field near the origin is of approximate solid-body (Rankine) 
form, if we subtract off this field and re-plot the homotopic sequence of patterns that emerges, shown in 
Figure |15[ a far richer and more instructive sequence of patterns is revealed, one that is far more relevant 
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Figure 13: Colored passive tracers advected by the velocity field z given in and depicted in the frame rotating with the 
vortex structure. As time evolves, the passive tracers stretch and mix forming large lobes at a finite distance from the initial 
location of the vortex structure. The separatices of the relative velocity field v — f?£^ are superimposed in black at various 
instants in time. 



for the understanding of the evolution of passive particle transport, as shown clearly in Figure [13| 

We finish by mentioning connections of this work in two other contexts. First, there is by now a 
growing body of work on calculating 'time-dependent separatrices' in developing flows that goes under the 
name of 'Lagrangian coherent structures' (LCS) [37JGEI- Certainly these tools are potentially useful for 
further elucidating the intermediate timescale dynamics associated with the evolution towards the Lamb- 
Oseen state, particularly for more complex initial patterns that perhaps start out as relative equilibria of 
the Euler equations. Second, if one regards, the vorticity field as a probability density function associated, 
for example, with the positions of initial system of point vortices undergoing a random walk, there are 
meaningful interpretations of the models used in this paper that have been discussed most recently, for 
example, in [SHI 001 SI] • While this interpretation has not been the main focus of our work, we do find it 
potentially ripe for future development. 
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